%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function generates a "smooth star" centered at xxc.
% The contour is parameterized as follows:
%
%   C(1,i) =                      x1 coordinate of node i
%   C(2,i) =        derivative of x1 coordinate of node i
%   C(3,i) = second derivative of x1 coordinate of node i
%   C(4,i) =                      x2 coordinate of node i
%   C(5,i) =        derivative of x2 coordinate of node i
%   C(6,i) = second derivative of x2 coordinate of node i
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function [C,curvelen,xx_int,xx_ext] = get_geometry_star(ntot,xxc,k,theta)

r        = 0.3;
%k        = 5;
tt       = linspace(0,2*pi*(1 - 1/ntot),ntot); 
C        = zeros(6,ntot );
% x(tt)
C(1,:)   =   1.5*cos(tt) + (r/2)*            cos((k+1)*tt + theta) + (r/2)* ...
    cos((k-1)*tt + theta);
% dx/dtt
C(2,:)   = - 1.5*sin(tt) - (r/2)*(k+1)*      sin((k+1)*tt + theta) - (r/2)* ...
    (k-1)*      sin((k-1)*tt + theta);
% d2x/dtt2
C(3,:)   = - 1.5*cos(tt) - (r/2)*(k+1)*(k+1)*cos((k+1)*tt + theta) - ...
    (r/2)*(k-1)*(k-1)*cos((k-1)*tt + theta);
% y(tt)
C(4,:)   =       sin(tt) + (r/2)*            sin((k+1)*tt + theta) - (r/2)* ...
    sin((k-1)*tt + theta);
% dy/dtt
C(5,:)   =       cos(tt) + (r/2)*(k+1)*      cos((k+1)*tt + theta) - (r/2)* ...
    (k-1)*      cos((k-1)*tt + theta);
% d2y/dtt2
C(6,:)   = -     sin(tt) - (r/2)*(k+1)*(k+1)*sin((k+1)*tt + theta) + ...
    (r/2)*(k-1)*(k-1)*sin((k-1)*tt + theta);

curvelen = 2*pi;

rmin = sqrt(min(C(1,:).^2 + C(4,:).^2));
rmax = sqrt(max(C(1,:).^2 + C(4,:).^2));

%%% Construct interior points.
ttint  = 2*pi*sort(rand(1,3));
xx_int = 0.5*rmin*[cos(ttint); sin(ttint)];

%%% Construct exterior points.
ttext  = 2*pi*sort(rand(1,5));
xx_ext = 1.5*rmax*[cos(ttext); sin(ttext)];

%%% Shift to center at xxc.
C([1,4],:) = C([1,4],:) + xxc*ones(1,size(C,2));
xx_int     = xx_int     + xxc*ones(1,size(xx_int,2));
xx_ext     = xx_ext     + xxc*ones(1,size(xx_ext,2));

% Rescale the derivatives in order to get the size of the element
h = curvelen / ntot;
C(2,:) = h* C(2,:);
C(5,:) = h* C(5,:);

return
